%% Saliency study - VAS analysis

cd('C:\Users\leu\OneDrive - UCL\Studies\Study - Stimulation Saliency\Registered Report\analysis\VAS analysis')

%% Get VAS data from individual MATLAB files

numfiles = 31; % 3 subject that only have 20 trials, run this section seperately for them and store data in matlab file
mydata = cell(1, numfiles);
N = cell(1, numfiles); %for saliency rating variable

highOddball_all_clean = [];
lowOddball_all_clean = [];
baseline_stim_peak= [];
highO_all= [];
lowO_all= [];

fstruct = dir('*_VAS*.mat');

for k = 1:numfiles
  myfilename = fstruct(k,1);
  mydata{k} = load(myfilename.name);
  baseline_stim_peak(:,k)= [mydata{k}.stim_temp];

   for m = 1:6      
      for f= 1: 4   
          trial_order =  string(mydata{k}.trial_order(m,f));

        switch  trial_order
            case 'highOddball'
              highOddball_all_clean(:,end+1) = mydata{k}.VAS_Results(:,(m-1)*f+f);
           
            case 'lowOddball'
             lowOddball_all_clean(:,end+1) = mydata{k}.VAS_Results(:,((m-1)*f+f));
            
        end
      end
   end
end

%% adding subjects with 20 and 24 trials together
cd('C:\Users\leu\OneDrive - UCL\Studies\Study - Stimulation Saliency\Registered Report\analysis\VAS analysis/not 24 trials/')
load('highOddball_not24.mat') %import data of subjects with not 24 trials from adequate folder
load('lowOddball_not24.mat')
load('baseline_stim_peak_not24.mat')

cd('C:\Users\leu\OneDrive - UCL\Studies\Study - Stimulation Saliency\Registered Report\analysis\VAS analysis')
highOddball_clean= [highOddball_all_clean, highOddball_not24]; 
lowOddball_clean= [lowOddball_all_clean, lowOddball_not24];
baseline_stim_peak_all = [baseline_stim_peak, baseline_stim_peak_not24];

% import cleaned variables in Letswave for FT analysis
% import MAT file: 
%                   Unit= amplitude, xUnit= time
%                   xStart=0, xStep = 1/ sampling rate = 0.01
%                   Dimension 1= x, Dimension 2= epoch

%% FIGURES
% plot low and high oddball condition together

time = (0.01:0.01:83)' ; %apostrophe transposes the vector to vertical

PeakSig=mean(highOddball_clean,2,'omitnan');
[pks,locs] = findpeaks(PeakSig(100:8300,1),time(100:8300,1),'MinPeakDistance',1.5);
pksH= pks(4:4:40);
locsH=locs(4:4:40);

PeakSig=mean(lowOddball_clean,2,'omitnan');
[pksL,locsL] = findpeaks(PeakSig(100:8300,1),time(100:8300,1),'MinPeakDistance',1.5); %find peak within certain distance to each other, -Variable to find negative peaks
locs_oddball = locsL(4:4:40);
pksLow = pksL(4:4:40);

plot(time, mean(highOddball_clean,2,'omitnan'),'LineWidth',3, 'Color','red', 'DisplayName','high oddball')
    xline([locsH],'--r','Color','red','LineWidth',1.1,'DisplayName','') %6s + 1.35s (as in Mulders and Colon, delay between 1.35 and 2s after stim)
hold on

plot(time, mean(lowOddball_clean,2),'LineWidth',3, 'Color',[0.00,0.45,0.74],'DisplayName','low oddball')
    xlabel('time [s]')
    ylabel('VAS rating [0-10]')
    xline(locs_oddball,'--r','LineWidth',1.1,'Color',[0.00,0.45,0.74],'DisplayName','')
    yline(5 ,'-.','Color','black')
    xlim([0,83])
    ylim([3,6])
    ax = gca;
    ax.FontSize = 25; 
    hA = get(gca);
    hA.XAxis.MinorTickValues = (2:2:83);
    hA.XAxis.MinorTick='on';
    box off

f=get(gca,'Children'); %get all legends into one vector
legend([f(23),f(12)],'high oddball','low oddball') %select only main curves, otherwise all horizontal and vertical lines are added as legends    

hold off

%% analyse oddball ratings

baseline_peak_mean = mean(baseline_stim_peak_all); % mean temp for baseline stimulation
baseline_peak_std = std(baseline_stim_peak_all); 

assessment_times = [7, 15, 23, 31, 39, 47, 55, 63, 71, 79];
meanCycleHighO = mean(locsH-assessment_times');
stdCycleHighO = std(locsH-assessment_times');
meanCycleLowO = mean(locs_oddball-assessment_times'); % compute the mean interval between maxima
stdCycleLowO = std(locs_oddball-assessment_times');

[h p ci stats] = ttest((locsH-assessment_times'), (locs_oddball-assessment_times')) %testing if timing of rating differs between conditions

avg_amp_VAS_highO=  mean(pksH); % from group level avg
std_amp_VAS_highO= std(pksH);
avg_amp_VAS_lowO= mean(pksLow);
std_amp_VAS_lowO= std(pksLow);

[h p ci stats] = ttest(pksH, pksLow) %testing if peaks at oddball are sign different from each other

%baseline
baseline_ampsHigh= mean(setdiff(pks, pksH,'stable')); %setdiff removes all rows that contain pksH from pks
baseline_ampsLow= mean(setdiff(pksL, pksLow,'stable'));
baseline_ampsHighstd= std(setdiff(pks, pksH,'stable')); 
baseline_ampsLowstd= std(setdiff(pksL, pksLow,'stable'));

%%
%test first vs third baseline cycle for potential differences
baseline_pksH= setdiff(pks, pksH,'stable');
baseline_pksL=setdiff(pksL, pksLow,'stable');

baseline_pksH_first= baseline_pksH(1:3:end);
baseline_pksH_third= baseline_pksH(3:3:end);
baseline_pksL_first= baseline_pksL(1:3:end);
baseline_pksL_third= baseline_pksL(3:3:end);

[h p ci stats] = ttest(baseline_pksH_first, baseline_pksH_third)
[h p ci stats] = ttest(baseline_pksL_first, baseline_pksL_third)

mean(baseline_pksH_first)
mean(baseline_pksH_third)
mean(baseline_pksL_first)
mean(baseline_pksL_third)

%% plot mean for each subject seperately
% %
% Oddball_pksH = [];
% Oddball_pksL = [];
% 
% for m = 1:31
%     [pksH,locsH] = findpeaks(highOddball_SUBmean(1:8300,m),time(1:8300,1),'MinPeakDistance',7);
%     figure
%    % findpeaks(highOddball_SUBmean(1:8300,m),time(1:8300,1),'MinPeakDistance',7);
%     [pksL,locsL] = findpeaks(lowOddball_SUBmean(100:8300,m),time(100:8300,1),'MinPeakDistance',1.5);
%     %findpeaks(lowOddball_SUBmean(100:8300,m),time(100:8300,1),'MinPeakDistance',1.5);
%     locs_oddball = locsL(4:4:end);
%     plot(time(100:8300,1), lowOddball_SUBmean(100:8300,m))
%         xline(locs_oddball,'--r')
% 
%     %Oddball_pksH(:,m)=pksH ;
%     %Oddball_pksL(:,m)=pksL(4:4:end);
% 
%     clear pksL
%     clear locsL
% end
% plot(time, highOddball_SUBmean)
% plot(time, lowOddball_SUBmean)

%% get rating peak for each participant within 1-2.6s after oddball stimulation peak, aggregate peaks
% Define the parameters
sampling_rate = 100; % in Hz
time_window = [1 2.6]; % in seconds
assessment_times = [7, 15, 23, 31, 39, 47, 55, 63, 71, 79]; % in seconds
col_num= width(highOddball_clean);

% Initialize an empty array to store the peak values
peak_valuesHigh = zeros(length(assessment_times), col_num);
peak_valuesLow = zeros(length(assessment_times), col_num);

%Loop over each participant
for col= 1:col_num
    % Loop over each assessment time
    for i = 1:length(assessment_times)
        % Define the start and end times for the window
        start_time = assessment_times(i) + time_window(1);
        end_time = assessment_times(i) + time_window(2);
        
        % Convert time to indices (sampling rate is 100 Hz, so 100 samples per second)
        start_index = round(start_time * sampling_rate);
        end_index = round(end_time * sampling_rate);
        
        % Extract the data within the time window
        window_dataHigh = highOddball_clean(start_index:end_index,col);
        window_dataLow = lowOddball_clean(start_index:end_index,col);

        % Find the peak value in the current window
        peak_valuesHigh(i,col) = max(window_dataHigh);
        peak_valuesLow(i,col) = max(window_dataLow);

    end
end

aggregateHigh = zeros(col_num,1);
aggregateLow = zeros(col_num,1);
for aggregate = 1: col_num
    aggregateHigh(aggregate,:)= (sum(peak_valuesHigh(:,aggregate))/10)';
    aggregateLow(aggregate,:)= (sum(peak_valuesLow(:,aggregate))/10)';
end

%export data for LMM in R
mean(aggregateHigh)
std(aggregateHigh)
mean(aggregateLow)
std(aggregateLow)


time = (0.01:0.01:83)' ; %apostrophe transposes the vector to vertical
plot(time, mean(highOddball_clean,2,'omitnan'),'LineWidth',3, 'Color','red')
 yline(mean(peak_valuesHigh,2))

plot(time, mean(lowOddball_clean,2,'omitnan'),'LineWidth',3, 'Color','blue')
 yline(mean(peak_valuesLow,2))


%% same extraction but for the peaks following baseline stimulation cycles

% Define the parameters
sampling_rate = 100; % in Hz
time_step = 2; % Step size for time intervals (2 seconds)
time_start = 1; % Starting time (1 second)
excluded_times = [7, 15, 23, 31, 39, 47, 55, 63, 71, 79]; % Times to exclude
time_window = [1 2.6]; % Window from 1 to 2.6 seconds after each time point
col_num= width(highOddball_clean);


% Generate the time points (1, 3, 5, 9, 11, 13, ...)
max_time = 80; % Maximum time for the assessment
all_times = time_start:time_step:max_time;
time_points = setdiff(all_times, excluded_times); % Remove the excluded times
num_columns = width(highOddball_clean); % Number of columns (datasets)

% Initialize an empty array to store the peak values
% Rows correspond to each time point, columns to each dataset
peak_valuesBase_high = zeros(length(time_points), num_columns);
peak_valuesBase_low = zeros(length(time_points), num_columns);

% Loop over each dataset (column) and each time point
for col = 1:num_columns
    for i = 1:length(time_points)
        % Define the start and end times for the window after the time point
        start_time = time_points(i) + time_window(1);
        end_time = time_points(i) + time_window(2);
        
        % Convert time to indices (sampling rate is 100 Hz, so 100 samples per second)
        start_index = round(start_time * sampling_rate);
        end_index = round(end_time * sampling_rate);
        
        % Extract the data within the time window for the current column
        window_dataBH = highOddball_clean(start_index:end_index, col);
        window_dataBL = lowOddball_clean(start_index:end_index, col);
        
        % Find the peak value in the current window and store it
        peak_valuesBase_high(i, col) = max(window_dataBH);
        peak_valuesBase_low(i, col) = max(window_dataBL);

    end
end

%aggregate all peaks (sum/num peaks), number peaks = 30
aggregateBase_high = zeros(col_num,1);
aggregateBase_low = zeros(col_num,1);
for aggregate = 1: col_num
    aggregateBase_high(aggregate,:)= (sum(peak_valuesBase_high(:,aggregate))/30)';
    aggregateBase_low(aggregate,:)= (sum(peak_valuesBase_low(:,aggregate))/30)';
end

%average all peaks 
aggregateBase_high = zeros(col_num,1);
aggregateBase_low = zeros(col_num,1);
for aggregate = 1: col_num
    aggregateBase_high(aggregate,:)= mean(peak_valuesBase_high(:,aggregate))';
    aggregateBase_low(aggregate,:)= mean(peak_valuesBase_low(:,aggregate))';
    aggregateBase_first(aggregate,:)= mean(peak_valuesBase_high(1:3:end,aggregate))'; %extract first cycle after oddball
    aggregateBase_third(aggregate,:)= mean(peak_valuesBase_high(3:3:end,aggregate))';%extract third cycle after oddball
end

%export data for LMM

mean(aggregateBase_high)
std(aggregateBase_high)
mean(aggregateBase_low)
std(aggregateBase_low)

%compare 1st vs 3rd cycle after oddball
mean(aggregateBase_first)
std(aggregateBase_first)
mean(aggregateBase_third)
std(aggregateBase_third)

[h p ci stats] = ttest(aggregateBase_first, aggregateBase_third, paired=true)


time = (0.01:0.01:83)' ; %apostrophe transposes the vector to vertical
plot(time, mean(highOddball_clean,2,'omitnan'),'LineWidth',3, 'Color','red')
 yline(mean(peak_valuesBase_high,2))

plot(time, mean(lowOddball_clean,2,'omitnan'),'LineWidth',3, 'Color','blue')
 yline(mean(peak_valuesBase_low,2))
